Hydroxide Diffusion in Functionalized Cylindrical Nanopores as Idealized Models of Anion Exchange Membrane Environments: An Ab Initio Molecular Dynamics Study

Anion exchange membranes (AEMs) have attracted significant interest for their applications in fuel cells and other electrochemical devices in recent years. Understanding water distributions and hydroxide transport mechanisms within AEMs is critical to improving their performance as concerns hydroxide conductivity. Recently, nanoconfined environments have been used to mimic AEM environments. Following this approach, we construct nanoconfined cylindrical pore structures using graphane nanotubes (GNs) functionalized with trimethylammonium cations as models of local AEM morphology. These structures were then used to investigate hydroxide transport using ab initio molecular dynamics (AIMD). The simulations showed that hydroxide transport is suppressed in these confined environments relative to the bulk solution although the mechanism is dominated by structural diffusion. One factor causing the suppressed hydroxide transport is the reduced proton transfer (PT) rates due to changes in hydroxide and water solvation patterns under confinement compared to bulk solution as well as strong interactions between hydroxide ions and the tethered cation groups.


■ INTRODUCTION
Anion exchange membrane (AEM) fuel cells (AEMFCs) rank among the cleanest low-cost electrochemical device technologies. The low cost is a result of the fact that AEMFCs operate in an alkaline environment, which obviates the need for precious metal catalysts 1−9 and gives AEMFCs an advantage over present proton exchange membrane fuel cells. Early versions of AEMFCs were used by NASA to power onboard systems during the Apollo and space-shuttle missions. 10 However, the need to enhance both hydroxide conductivity and chemical stability in AEMFCs remains an impediment to their widespread adoption. Designing AEMs with high hydroxide conductivity and chemical stability requires a detailed understanding of the hydroxide transport mechanism in geometries reflective of AEM environments.
Molecular dynamics (MD) is a powerful tool for studying the atomistic mechanisms that govern hydroxide transport mechanisms in both bulk and nanoconfined environments. Hydroxide diffuses structurally (also known as the Grotthuss mechanism) via individual sequential proton transfer (PT) events initiated by presolvation of the hydroxide ion. 11−17 Previous MD studies of hydroxide transport in AEM systems employed reactive force fields to achieve chemical bondbreaking and forming in PT. 5,18−22 A multistate empirical valence bond (MS-EVB) study of a poly(vinyl benzyltrimethylammonium) (PVBTMA) system indicated that the hydroxide transported within the first solvation shell of trimethylammonium (TMA) cations with a greater contribution from vehicular mechanism. The regularly spaced TMA cations observed in this system form overlapping first solvation shell regions, allowing hydroxide to transport from the solvation shell of one cation to the other. In comparison, a series of ReaxFF simulations carried out in functionalized poly(phenylene oxide) (PPO) systems suggested that the structural diffusion mechanism is dominant and, in particular, is critical for hydroxide ions to pass through bottlenecks of water channels where vehicular diffusion is associated with high energy barriers needed to shed solvating waters upon entering these bottleneck regions. 18−21,23−26 The dominance of structural diffusion is also supported by various experimental measurements 27 as well as overly low vehicular diffusion coefficients of hydroxide from nonreactive force field MD simulations. 4,18−21,28−32 Ab initio MD 33−35 (AIMD) simulations have proved enormously successful in predicting hydroxide diffusion constants, mechanisms, and kinetics in aqueous solution and in revealing both the presolvation concept and the dynamical hyper-coordination picture of hydroxide transport. 13,14,17,36−40 In an AIMD simulation, nuclear motion is generated from Newton's second law of motion using forces generated "on the fly" from electronic structure calculations. The high computational overhead and O(N 3 ) system size scaling of AIMD precludes modeling AEMs at the mesoscale. On the other hand, nanoconfined environments have been utilized in both experimental and theoretical studies to model polymer architectures and nanoporous materials for use in electrochemical devices. 23−26,41−51 Materials with cylindrical pores, in particular, have attracted interest for their ability to enhance diffusion in proton-conducting liquids. 41,44,45,50,52−54 In addition, polymers such as PPO can support morphologies containing tight pores for different linker lengths to the TMA cations. 55 −57 In order to capture the effect of the AEM nanoconfined environment on hydroxide diffusion mechanisms, idealized geometries such as graphane bilayers and carbon nanotubes to which cationic functional groups are attached have proved successful. 23−26,49,51 Our previous work employed graphane bilayer structures as an idealized mimic of a lamellar AEM morphology. 23−26,49,51 In this work, we consider an idealized mimic of the type of nanopore that might arise in bicontinuous morphologies by employing graphane nanotube (GN) structures functionalized by cations attached to the inside of the tube via short linkers. AIMD simulations are performed in order to probe the influence of the cylindrical confinement and cation placement. The simulations indicate that although dominated by structural diffusion, the hydroxide transport is suppressed in these confined systems compared to bulk solution. We attribute the suppressed hydroxide transport to the decreased PT rate due to the high barriers associated with the presolvation under nanoconfinement, as well as strong hydroxide−cation interactions.

■ METHODS
Snapshots depicting our GN models are provided in Figure 1, and system parameters are listed in Table 1. Within this model, zigzag GNs GN(n,0) mimic the hydrocarbon backbone in a local cylindrical pore structure in an AEM. The GNs are aligned along thez axis with periodic boundary conditions applied in this direction to form channels of   64 and molecular sieves. 65 AIMD simulations of these confined systems were performed with the CPMD code. 66 The electron structure is described at the level of density functional theory (DFT) 67−69 using the BLYP generalized gradient functional. 70,71 Despite the approximations inherent in the use of the BLYP functional for aqueous systems, our choice of this functional for this study is based on several considerations. First, previous studies with it demonstrated its ability to produce correct hydroxide transport mechanisms in bulk aqueous solutions 14,36−38,72−74 and to yield new and testable predictions of hydroxide transport under nanoconfinement that have been confirmed experimentally. 26,75 Moreover, when used with a converged basis set, its description of the structural and dynamical properties of bulk water are reasonable. 76−78 In order to connect the present study with our previous work, we, therefore, opted to employ the BLYP functional again here; however, in recognition of the need for a large basis set, we used a cutoff value of 80 Ry for the plane wave basis set. We are, however, aware of the caveats associated with the use of this functional and will address improvements for future studies in the Conclusions section. In this study, dispersion interactions were included via the dispersion-corrected atomic core pseudopotentials (DCACP) approach. 79,80 The Car−Parrinello algorithm 81 is used with a fictitious electron mass of 600 a.u. and a time step of 4 a.u. All GN carbon and hydrogen coordinates were kept fixed throughout the simulations. The deuterium mass was used for all hydrogen atoms to reduce the importance of nuclear quantum effects. After wave function minimization, each system was heated to approximately 300 K and then equilibrated in the NVT ensemble with massive Nose−Hoover chain thermostats 82 (12 ps for GN(16,0), 13 ps for GN(17,0), and 10 ps for GN(19,0)). Following equilibration, production runs were performed in the NVE ensemble (87 ps for GN (16,0), 109 ps for GN (17,0), and 73 ps for GN (19,0)). These simulation settings are consistent with our previous studies in systems confined between graphane bilayers. 23−26,49,51 ■ RESULTS AND DISCUSSION Diffusion Coefficients of Hydroxide. The diffusion coefficients of hydroxide ions serve as the microscopic indicator of hydroxide conductivity. Structural diffusion causes the identity of the hydroxide to change with each PT reaction. Thus, in order to identify hydroxide ions in each AIMD step, each solution-phase hydrogen atom is assigned to its closest oxygen atom. Oxygen atoms with only one assigned hydrogen are recognized as hydroxide, for which the oxygen atom is denoted O* and the corresponding hydrogen as H*. Oxygen atoms with two assigned hydrogen atoms are water oxygen atoms O W with hydrogen atoms denoted H W .
Hydroxide diffusion coefficients are calculated from the mean square displacement (MSD) of the O* coordinates, including all changes of O* identity due to PT (see Figure S1 of the Supporting Information). The results, together with water O W diffusion coefficients (corrected for finite-size effects, 77,86 ) are given in Table 2. The GN(19,0) system has the highest hydroxide and water diffusion coefficients, while the GN(17,0) system has the lowest hydroxide and water diffusion coefficients. Note that even the highest hydroxide and water diffusion coefficients in GN (19,0) are still smaller than the previously reported values in bulk aqueous solution obtained from experiments 14,38,83 and AIMD simulations with the same functional. 74,77 The GN(16,0) system with a lower water content (λ = 10) has a smaller water diffusion coefficient than the GN(19,0) system, but the hydroxide diffusion coefficients are close in these two systems. This retention of hydroxide diffusivity under decreased water diffusivity suggests a hydroxide transport mechanism dominated by structural diffusion in GN(16,0). The further reduced hydroxide diffusivity and water diffusivity in GN(17,0) results from the high solution-phase density in this system.
Structure of the Solution Phase under Nanoconfinement. Understanding the solution-phase structure under nanoconfinement is fundamental for further analysis of the hydroxide transport mechanism. The geometrical analysis will make use of standard cylindrical coordinates (ρ, ϕ, z), where ρ is the radial coordinate, ϕ is the azimuthal angle, and z is the coordinate along the length of the GN. Distributions of oxygen and hydrogen atoms in Figure 2 are calculated using 1 (below), which prescribes the probability density as a function of the cylindrical radius ρ (ρ = 0 corresponds to the GN axis): The bracket represents the ensemble average over all configurations. The summation is taken over all the corresponding atoms. The bin width for the statistics is δρ = 0.05 Å. z GN is the periodic length of the GNs. n is the numerical density of the corresponding species in the solution phase and calculated with the available volume inside the GN (V conf ). Complementary analysis of oxygen distributions along the ϕ and z coordinates is provided in Section S2 of the Supporting Information. In Figure 2a, the oxygen distribution peaks at ρ > 2 Å in all three systems suggest that the solution forms cylindrical layers near the GN wall. The radii of these cylindrical layers reflected in the peak positions increase with GN radius. Similar cylindrical layer structures were also observed for pure water 87 and aqueous acidic solutions 88,89 confined in carbon nanotubes, and they also closely resemble the layered structures seen in our confined systems between graphane bilayers. 24,26,49 However, since we explicitly included cations in our nanoconfined system, the solution phase is mostly located in the region of − 90°< ϕ < 90°. In other words, the cylindrical layers are incomplete. Oxygen distributions are also observed at small ρ values (near the GN axes) in GN(17,0) and GN(19,0). As illustrated in Figure 1d,1g, these oxygen densities do not correspond to cylindrical structures. The main peaks of the hydrogen density distribution reside slightly inside the corresponding oxygen cylindrical layers. The small peaks or shoulders in the hydrogen distribution outside the oxygen cylindrical shell (ρ =3.5 Å for GN(16,0), 3.9 Å for GN(17,0), and 4.7 Å for GN(19,0)) suggest dangling hydrogen atoms that form no hydrogen bonds (HBs) but instead point toward the GN wall.  91,94 In the next section, we will show that this radially aligned orientation and the absence of donating HBs (dangling hydrogen) of this O*−H* bond ultimately affect the hydroxide PT mechanism.
Detailed Mechanism of Hydroxide PT Reaction under Confinement. PT is the fundamental step in the structural diffusion process of hydroxide ions, and therefore, the PT mechanism affects the hydroxide transport rate. In Figure 4a, the free energy of PT reaction is calculated as δ min is the minimum over all δ values among the HBs accepted by a hydroxide anion. The HB with δ min is designated as the most active HB in the system, meaning that PT through this HB is most likely. The water oxygen in this HB is denoted by O W (1 * ) . A small δ min (< 0.1 Å) corresponds to a proton that is nearly equally shared between the two oxygen atoms, i.e., the PT transition state. A large δ min (> 0.5 Å) value indicates that the proton is attached to one of the oxygen atoms and can be considered as a PT resting state. The free energy profiles are symmetrized about δ min = 0, and minima are shifted to zero. All resulting free energy profiles are double-well shaped with minima at δ min = ± 0.46 Å and barriers at δ min = 0. GN(19,0) and GN(16,0) and a continuous (history dependent) population correlation function as The function C i (t) gives the probability that an oxygen that was O* at t = 0 becomes O* again at time t independent of any changes in identity between 0 and t. The function C c (t) gives the probability that O* retains its identity for a time t. One process contributing significantly to C c (t) is proton rattling, i.e., the proton shuttles back and forth along the same HB, which contributes nothing to overall hydroxide transport. The long timescale behavior of the hydroxide identity decay (slow processes) are contained in the intermittent population correlation function C i (t). As shown in Figure 4b, for both the continuous and the intermittent population correlation functions, GN(17,0) has the slowest hydroxide identity decay, while GN(19,0) has the fastest. This ordering of the hydroxide identity decay rate (PT rate) is consistent with the ordering of the hydroxide and water diffusion coefficients in Table 2.
The PT rate is closely related to the underlying molecular mechanism. In the dilute bulk aqueous solution, PT of hydroxide follows the dynamic hyper-coordination mechanism. 13,14,17,38 Hydroxide at the PT resting state is preferentially hyper-coordinated, i.e., accepting four HBs from solvating water molecules and occasionally donating one. During a PT reaction, the hydroxide loses an accepting HB and strengthens the donating HB. Thus the hydroxide forms a water-like tetrahedral solvation structure with three accepting HBs and one donating HB and is, thus, properly "pre-solvated" for a subsequent PT reaction. However, PT via the "PT-active complex" (a tetrahedrally coordinated hydroxide and proton-donating water) at the transition state is not the only feasible pathway, and, in fact, changes in the local environment dictate which pathway dominates the PT process. Recent work reported simulations of aqueous NaOH solutions using a neural network potential. 95 The "most important", or most active, presolvation structure, defined as the structure with the highest PT rate and the lowest free energy barrier, changes with NaOH concentration. These authors attributed this observation to changes in the solvation configurational distributions of hydroxide and water with increasing NaOH concentrations. Accordingly, the presolvation process, which requires a change from the most probable solvation configuration to the most active configuration, was achieved by changes in accepting HBs, donating HBs, or coordinating Na + cations of the hydroxide or the proton-donating first solvation shell water. Furthermore, we noticed that in all the most active presolvated PT complexes, hydroxide and proton-donating water were similarly solvated and produced symmetric PT free energy profiles.
Here we will show that similar to the NaOH concentration dependence, the confined environment of the GN also modulates the solvation configurational distributions of hydroxide and water, prompting different hydroxide PT pathways and the properly presolvated PT complexes. Solvation configurational probability distributions (Table 3) were calculated as the probability of finding a hydroxide or water molecule that donates n HBs and accepts m HBs (noted as nDmA  Table 4).   (16,0) is associated with the breaking of an accepted HB of the hydroxide, in agreement with the dynamic hyper-coordination mechanism in the dilute bulk solution. However, the occasional donating HB of hydroxide in bulk aqueous solution is almost absent in the GN(16,0) system as indicated by the vanishing probabilities of 1DmA hydroxide configurations at both the PT resting state and the PT transition state, and the presolvated hydroxide configuration for PT is the trigonal pyramidal 0D3A (71.1%) instead of the full tetrahedral 1D3A structure (0.6%). Correspondingly, the most active first solvation shell water at the PT transition state is dominated by a 1D2A pattern (71.3%), the configuration also involving a dangling hydrogen. This substantial population of 1D2A at the PT transition state is in stark contrast with the low population (17.0%) of 1D2A for all water molecules in the solution phase, which can be attributed to the asymmetry between water donor and acceptor sites. 96 Therefore, in addition to breaking an accepting HB of the hydroxide, presolvating the PT complex should also include an extra step in which the protondonating water adopts the 1D2A configuration. This associated energy cost of this extra step reduces the PT rate and provides an explanation for the suppressed hydroxide diffusivity within GN confinement. Similar suppression of hydroxide transport is observed in a previous study of hydroxide migration in a confined water monolayer between graphene sheets. 97 In that study, the proton-donating water must also adopt the low population 1D2a pattern to be presolvated for PT to hydroxide.
As a water molecule assumes the 1D2A solvation configuration, we also observe reorientation of its O−H bond to be aligned along the radial direction and to point to the GN wall. This reorientation is captured in the probability distribution of the O W −H W − radial angle θ at different stages of the PT reaction in Figure 5.  GN(17,0). Nevertheless, the high solution density in GN(17,0) creates another PT pathway. The very high population (90.0%) of the 0D4A pattern and low (7.2%) population of the 0D3A hydroxide solvation structure indicate that it is highly energetically unfavorable to break an accepting HB from the hyper-coordinated hydroxide solvation complex. Instead, PT may involve a hypercoordinated 0D4A hydroxide pattern and a similarly solvated "hyper-coordinated" 1D3A water as evidenced by the high (57.2%) 0D4A population for hydroxide and the second-  highest (37.6%) 1D3A population for the most active water at the PT transition state. The hyper-coordinated water, whose oxygen atom has also been referred to as an overcoordinated oxygen (OCO), 96 has one hydrogen participating in the HB donated to the hydroxide while the other hydrogen dangles. Similar most important or most active PT complexes involving this hyper-coordinated water were also reported in the study of aqueous NaOH solution systems. 95 This "unconventional" solvation structure of the PT complex involving a hyper-coordinated hydroxide and a hypercoordinated water provides an explanation for the higher PT reaction free energy barrier, lower PT rate, and further reduced hydroxide diffusivity of the GN(17,0) system. A detailed analysis of the hyper-coordinated water is provided in Section S4 of the Supporting Information.
The dominant solvation configurations in GN (19,0) at the PT transition state are the 0D3A hydroxide pattern (49.9%) and the 1D2A complex for the most active water (51.4%), which are the same configurations as in GN(16,0). The weak HB donated by hydroxide is present in GN(19,0) and participates in PT reactions as indicated by the higher 1DmA hydroxide populations at the PT transition state in GN (19,0) than in GN(16,0) and GN (17,0). The corresponding 2DmA pattern for the most active water (without a dangling hydrogen) also has higher populations at the PT transition state than the other two systems. The orientational preference of the O W  Figure 5c. Overall, with the higher water content and larger GN radius, hydroxide PT in GN(19,0) exhibits features of the dynamic hypercoordination mechanism and more closely resembles dilute bulk aqueous solution behavior than is seen in GN(16,0) and GN(17,0).
Hydroxide Diffusion Characteristics in the GN Confined System. In addition to the fundamental PT mechanism, other factors, including competition with vehicular diffusion and interaction with the cations, also affect the hydroxide transport. Here we provide further analysis of these factors.
To quantify the relative contribution from structural diffusion and vehicular diffusion, MSD is decomposed into discrete and continuous components according to 5,98 The bracket represents the ensemble average over all configurations. MSD d and MSD c are the discrete and continuous components corresponding to structural and vehicular diffusion, respectively. The function ⟨Δr d (t) · Δr c (t)⟩ is the correlation between these two components. The discrete displacement Δr d (t) and the continuous displacement Δr c (t), for the period from time 0 to time t, are calculated by accumulating the corresponding elementary displacements δr d (t, δt) and δr c (t, δt). If PT occurs during the period from t to t + δt, the hydroxide displacement during this time window δt is assigned to the elementary discrete displacement δr d (t, δt) and 0 is assigned to the elementary continuous component δr c (t, δt). Otherwise, the hydroxide displacement is assigned to the elementary continuous component δr c (t, δt) and 0 is assigned to the elementary discrete component δr d (t, δt). The time window δt is selected to be approximately 97 fs, which is close to the 100 fs used in a previous study. 5 We noticed that the absolute amplitudes of the decomposed MSD curves are affected by the length of the time window, but the relative ordering of the contributions (i.e., discrete > continuous) remains the same (see Section S5 of the Supporting Information). As shown in Figure 6, the discrete components of all three systems contribute significantly more than the continuous components, indicating that structural diffusion dominates the hydroxide diffusion mechanism in all three systems. The dominance of structural diffusion in the hydroxide transport process explains why hydroxide diffusivity is maintained at lower water content and lower water diffusivity in GN (16,0). The anticorrelation between discrete and continuous components suggests a competition between structural and vehicular diffusion, which originates from the presence of The Journal of Physical Chemistry C pubs.acs.org/JPCC Article cations. When the hydroxide hops away from the cation via PT, vehicular diffusion can drag it back toward the cation, and when the hydroxide diffuses away from the cation via the vehicular mechanism, the structural diffusion can bring it back via PT. Similar anticorrelation between the discrete and continuous components was also reported in the MS-EVB simulation of PVBTMA AEM, 5 while almost zero correlation was observed for an excess proton in aqueous solution confined within a carbon nanotube with no counter ions. 98 The strong cation−hydroxide interaction can be illustrated by plotting the z-coordinates of the hydroxide O* and cation nitrogen as functions of time (see Figure 7). The figure shows that the hydroxides remain close to one cation within our simulation. One exception is the hydroxide of the red trajectory in GN (16,0) in Figure 7a, which undergoes two successful PT events (∼7 and ∼38 ps) and hops from one cation to the other. The other exception is the hydroxide in the red trajectory of GN(19,0) (Figure 7c), which mainly resides within the region between two cations.
Further analysis of the hydroxide−cation interaction is performed by dividing the first solvation shells of the TMA cations into two regions�the overlapping region of two first solvation shells of two neighboring cations and the nonoverlapping region of the first solvation shells belonging to just one cation. An N−O distance cutoff of 5.95 Å, corresponding to the first minimum of N−O RDFs (see Figure S7 of the Supporting Information), defines the first solvation shell of TMA cations. The spatial distribution of the overlapping and the nonoverlapping regions of the first solvation shells of cations is illustrated in Figure 8. The blue "voids" centered on ϕ = ± 180°are occupied by the cations. The red strips correspond to the overlapping region, and the white area corresponds to the nonoverlapping region. In GN(16,0) and GN (17,0), the space occupied by the solution under nanoconfinement only consists of overlapping and nonoverlapping regions. In other words, most of the water and hydroxide species are in the first solvation shell of at least one cation in these two systems. The second solvation shells (and beyond) of TMA cations only appear in GN(19,0) as the blue area centered on ϕ = 0°on the outer cylindrical layer (Figure 8c).
In order to make explicit the roles played by the overlapping and nonoverlapping regions of the first solvation shell of TMA cations in hydroxide transport, we calculated population distributions of the hydroxide O* (as well as all oxygen atoms O) in these two regions and listed them in Table 5. In GN(16,0) and GN(17,0), 99% water and hydroxide molecules exist in the first solvation shells of the cations, consistent with the previous analysis of spatial distributions of the overlapping and nonoverlapping regions. The hydroxide defect can visit different oxygen atoms over time under structural diffusion. If we assume that the hydroxide identity obeys a random distribution among all oxygen atoms, the distribution of hydroxide in the overlapping and nonoverlapping regions would be the same as that of all oxygen atoms. Instead, in both systems, the hydroxide O* has a higher population in the nonoverlapping first solvation shell regions than all oxygen atoms. In other words, the overlapping regions appear as barriers to hydroxide transport, consistent with the limited hydroxide transport from the vicinity of one TMA cation to that of the other in GN(16,0) and GN (17,0).
The GN(19,0) system has a larger GN radius and a higher water content, and consequently, the confined space within GN (19,0) accommodates an area beyond the first solvation shells of any cations. Compared to GN(16,0) and GN(17,0), the hydroxide O* in GN(19,0) has a decreased population in the nonoverlapping regions and an increased population in the overlapping regions, approaching the distribution of water in these two regions. Therefore, hydroxide diffuses more freely between these two regions in the GN(19,0) system. On the other hand, the total population of hydroxide O* in the first solvation shell of at least one cation is around 71%, while this population of all oxygen atoms is 61%, indicating that the hydroxide dissociation from the cations is still incomplete.

■ CONCLUSIONS
We constructed nanoconfined aqueous alkaline solution systems within graphene nanotubes functionalized with TMA cations as idealized models of cylindrical pores in AEMs and simulated these systems using AIMD to study how the hydroxide transport is affected by the confined  Correspondingly, water molecules in the first solvation shell water of hydroxide donating HBs to the hydroxide O* atoms also possess a dangling hydrogen oriented toward the GN wall during the PT, producing a symmetrically solvated PT complex. This unfavorable solvation configuration of water (donating just one HB) at the PT transition state reduces the PT rate and serves as one factor causing suppression of hydroxide diffusion in these confined systems. While this is the prominent PT pathway in GN(16,0), in GN(19,0), which has a higher water content, PT to the hydroxide ions more closely resembles that in bulk solution. In GN(17,0), the hyper-coordinated hydroxide and the symmetrically solvated hyper-coordinated water also participate in PT under the high solution density, further reducing the PT rate and hydroxide diffusivity in that system. The other factor impeding hydroxide transport is the strong attraction between hydroxide ions and the TMA cations. The hydroxide is preferentially associated with one cation instead of shared between two cations. And the hydroxide does not fully dissociate from the cation even at higher water content in GN (19,0).
Because of the dominance of discrete, structural diffusion over continuous, vehicular diffusion, as indicated in Figure 6, it is natural to ask about the role of nuclear quantum effects. The importance of nuclear quantum effects in bulk hydroxide transport was explored using ab initio path integral calculations in ref 36. There, it was seen that while the quantitative picture changes when nuclear quantum effects are included, the basic mechanistic picture does not. While it is difficult to speculate about how these effects would influence the quantitative results obtained here, it is not expected that the comparative picture of the hydroxide transport characteristics would change significantly by their inclusion. Similarly, because of the relative importance of structural diffusion, which is driven largely by local effects such as dynamical solvation patterns and the hydrogen-bond rearrangement, we would not expect the relatively small GN lengths to influence the mechanistic picture presented much Figure 8. Spatial distributions of the first solvation shells of cations calculated as the average numbers of cation first solvation shells in which an oxygen atom resides (i.e., average numbers of cations that an oxygen atom is contacting with). The angle ϕ is the azimuthal angle with cations attached to ϕ ≈ 0°, and z is the axial coordinate. In GN(19,0), a distance cutoff of r = 2.175 Å from the minimum of oxygen distributions along radius in Figure 2a is used to distinguish the outer cylindrical layer and the inner near-axis area. Red regions correspond to the overlapping areas between the first solvation shells of two neighboring cations. Gray regions correspond to the nonoverlapping regions of the first solvation shell of only one cation. The blue regions at ϕ = ± 180°are the spaces occupied by cations inside the confinement, while the blue region in (c) at ϕ = 0°corresponds to the solution phase beyond the first solvation shells of any cations. The Journal of Physical Chemistry C pubs.acs.org/JPCC Article as they are not seen to be significant in bulk aqueous solution. 14 With the understanding of hydroxide transport in the present nanoconfined model systems, we hope to provide insights that are helpful in the design of AEMs with higher hydroxide conductivity. Since the hydroxide diffusivity is suppressed in these confined systems, avoiding the formation of water channel bottlenecks is desired by minimizing the channel size distributions at small radii. If narrow water channels are inevitable, sizes down to our GN(16,0) (∼6.2 Å) are acceptable as the structural diffusion can help retain the hydroxide diffusivity. Changing the cation chemistry might also help enhance the structural diffusion component of the hydroxide transport, which we will investigate in future work. It is also desired to have complete dissociation between hydroxide and cations, which may be achieved by engineering linkers and polymer backbones so that the cations are buried inside the hydrophobic phases and separated from the hydrophilic phases.
In this study, the electronic structure calculations employed in the AIMD simulations involved significant approximations inherent in the BLYP functional. As Table 2 shows, these approximations affect the diffusivity of water at 300 K, causing the ratio of hydroxide to water diffusion constants to be too high. For the system sizes considered here and in our previous work, 23−26,49,51,75 employing a higher-level functional such as the SCAN meta-GGA functional, 99 which has been shown to soften the water structure compared to BLYP 100 and offer some improvements in the details of the solvation patterns of aqueous hydroxide, 101 would be worth exploring.
As an alternative to meta-GGA functionals, hybrid functionals could be employed. However, for the system sizes employed in this study, a time-saving feature, such as the combination of adaptively compressed exchange operator schemes and multiple time-stepping, 102,103 would be needed. Ultimately, extending the system size in order to examine more sophisticated and representative models will be desirable, at which point, AIMD simulations will become infeasible, and more efficient approaches such as reactive neural network potentials 104 become an attractive option that can retain the accuracy of AIMD if they are trained appropriately. ■ ASSOCIATED CONTENT
MSD plots (Section S1); oxygen distributions along the azimuthal angle ϕ and axial coordinate z (Section S2); RDFs and angle distributions for defining HB criteria and cation solvation shells (Section S3); detailed analysis of the hyper-coordinated water (Section S4); and decomposition of the MSD into discrete and continuous compositions with a short time window (Section S5) (PDF)